home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / MSDOS / (m)aaj / PLANETS.PAS < prev    next >
Pascal/Delphi Source File  |  1985-10-18  |  31KB  |  1,076 lines
  1. PROGRAM PLANETS;
  2. {
  3.  
  4.                                Version 1.0
  5.  
  6.      This Program computes information relating to the position, distance,
  7.   magnitude, etc for the major planets on a specified date and time.
  8.   Refer to PRACTICAL ASTRONOMY WITH YOUR CALCULATOR by Peter Duffett-Smith
  9.   for the calculation methods.
  10.  
  11.      This program is placed in the public domain "as is".
  12.  
  13.      Comments may be sent to the author:
  14.  
  15.              Larry Puhl
  16.              6 Plum Court
  17.              Sleepy Hollow, Ill. 60118
  18. }
  19.  
  20. CONST
  21.    W = 12; {Width of real printouts}
  22.    D = 6;  {Decinal places in real printouts}
  23.    Number_of_Planets = 9;
  24.    CenterX = 300;
  25.    CenterY = 92;
  26.    Screen_Aspect = 2.7;
  27.    Scale_for_inner_planets = 60;
  28.    Scale_for_outer_planets = 2;
  29.  
  30.    Display_Degrees_in_decimal : boolean = false;
  31.    Daylight_Savings_Correction_Enabled : boolean = true;
  32.    Display_Hour_Angle_in_Degrees : boolean = false;
  33.    Display_RA_in_Degrees : boolean = false;
  34.    Zone_correction : integer = 6;
  35.    Longitude : real = 88.3;
  36.    Latitude : real = 43.1;
  37.  
  38.  
  39. TYPE
  40.    orbit =
  41.       RECORD
  42.          name : STRING[10];
  43.          Epoch : integer;
  44.          Period : real;
  45.          Long_at_Epoch : real;
  46.          Long_of_Per : real;
  47.          Ecc : real;
  48.          Axis : real;
  49.          Inc : real;
  50.          Long_of_Ascend_Node : real;
  51.          Size : real;
  52.          Brightness : real;
  53.       END;
  54.    MSDOS_REGISTERS =
  55.       RECORD
  56.          AX,BX,CX,DX,BP,SI,DI,DS,ES,FLAGS: Integer;
  57.       END;
  58.  
  59. VAR
  60.    Bios_Display_Mode : Integer Absolute $0000 : $0465;
  61.    MSDOS_time : MSDOS_REGISTERS;
  62.    ch : char;
  63.    zero : real; {used to print 0}
  64.    planets : ARRAY [1 .. 10] OF orbit;
  65.    num : integer;
  66.    day,month,year,day_of_year : integer;
  67.    JD : real;
  68.  
  69.  
  70.    planet_long_of_ascend_node : real;
  71.    planet_mean_anomaly : real;
  72.    planet_heliocentric_long : real;
  73.    planet_true_anomaly : real;
  74.    planet_radius_vector : real;
  75.    planet_helio_ecliptic_lat : real;
  76.  
  77.    earth_long_of_ascend_node : real;
  78.    earth_mean_anomaly : real;
  79.    earth_heliocentric_long : real;
  80.    earth_true_anomaly : real;
  81.    earth_radius_vector : real;
  82.  
  83.    GMT : real; {Greenwich mean time}
  84.    GST : real; {Greenwich Siderial time}
  85.    LCT : real; {Local Civil time}
  86.    LST : real; {Local Siderial time}
  87.  
  88.    Hour : integer;
  89.    Minute : integer;
  90.    T,B : real;
  91.  
  92.    projected_heliocentric_long : real;
  93.    projected_radius_vector : real;
  94.  
  95.    Geocentric_ecliptic_Long : real;
  96.    Geocentric_latitude : real;
  97.  
  98.    RA : real;{equatoriprojected_radius_vector : real;
  99.  
  100.    Geocentric_ecliptic_Long : real;
  101.    Geocentric_latitude : real;
  102.  
  103.    RA : real;{equatorial coordinates}
  104.    DEC : real;
  105.  
  106.    Phase : real;
  107.    Distance_from_Earth : real;
  108.    Diameter : real;
  109.    Magnitude : real;
  110. al coordinates}
  111.    DEC : real;
  112.  
  113.    Phase : real;
  114.    Distance_from_Earth : real;
  115.    Diameter : real;
  116.    Magnitude : real;
  117.  
  118.    Hour_Angle : real;
  119.    Hour_Angle_in_Degrees : real;
  120.    Azimuth : real;
  121.    Altitude : real;
  122.  
  123.    LST_Rise : real;
  124.    LST_Set  : real;
  125.    GST_Rise : real;
  126.    GST_Set  : real;
  127.    GMT_Rise : real;
  128.    GMT_Set  : real;
  129.    LCT_Rise : real;
  130.    LCT_Set  : real;
  131.    Azimuth_Set : real;
  132.    Azimuth_Rise : real;
  133.  
  134.    y : real; {intermediate variable}
  135.    x : real; {intermediate variable}
  136.    z : real; {intermediate variable}
  137.    zmax : real; {intermediate variable}
  138.    time : real;
  139.    Daylight_savings_start : real;
  140.    Daylight_savings_end  : real;
  141.  
  142.  
  143. {       Trig Functions     }
  144.  
  145.  
  146.    FUNCTION DEGREE (X:real) : real;
  147.       BEGIN
  148.          DEGREE := X*57.295779513
  149.       END; {DEGREE}
  150.  
  151.    FUNCTION RADIAN (X:real) : real;
  152.       BEGIN
  153.          RADIAN := X/57.295779513
  154.       END; {RADIAN}
  155.  
  156.    FUNCTION SIND (X:real) : real;
  157.       BEGIN
  158.          SIND := SIN(X/57.295779513)
  159.       END; {SIND}
  160.  
  161.    FUNCTION COSD (X:real) : real;
  162.       BEGIN
  163.          COSD := COS(X/57.295779513)
  164.       END; {COSD}
  165.  
  166.    FUNCTION TAN (X:real) : real;
  167.       BEGIN
  168.          TAN := SIN(X)/COS(X)
  169.       END; {TAN}
  170.  
  171.    FUNCTION TAND (X:real) : real;
  172.       BEGIN
  173.          TAND := TAN(X/57.295779513)
  174.       END; {TAND}
  175.  
  176.    FUNCTION ARCTAND (X:real) : real;
  177.       BEGIN
  178.          ARCTAND := 57.295779513*ARCTAN(X)
  179.       END; {ARCTAND}
  180.  
  181.    FUNCTION ARCTANYX(Y,X : real) : real;
  182.       VAR
  183.  EGIN
  184.          ARCTAND := 57.295779513*ARCTAN(X)
  185.       END; {ARCTAND}
  186.  
  187.    FUNCTION ARCTANYX(Y,X : real) : real;
  188.       VAR
  189.          z,zmax : real;
  190.       BEGIN
  191.          z :=ArcTanD(y/x);
  192.          IF (y>0) AND (x>0) THEN zmax := 90;
  193.          IF (y>0) AND (x<0) THEN zmax := 180;
  194.          IF (y<0) AND (x<0) THEN zmax := 270;
  195.          IF (y<0) AND (x>0) THEN zmax := 360;
  196.          WHILE z>zmax DO z := z -180;
  197.          WHILE z<(zmax-90) DO z := Z + 180;
  198.          ARCTANYX := z;
  199.       END; {ARCTANYX}
  200.  
  201.    FUNCTION ARCSIN (X:real) : real;
  202.       BEGIN
  203.          IF X * X >= 1 THEN ARCSIN := 1.570796327 ELSE
  204.          ARCSIN := ARCTAN(X/SQRT(1-X*X))
  205.       END; {ARCSIN}
  206.  
  207.    FUNCTION ARCSIND (X:real) : real;
  208.       BEGIN
  209.          ARCSIND := 57.295779513 * ARCSIN(X)
  210.       END; {ARCSIND}
  211.  
  212.    FUNCTION ARCCOS (X:real) : real;
  213.       BEGIN
  214.          IF X * X >= 1 THEN ARCCOS := 1 ELSE
  215.          ARCCOS := 1.570796327-ARCTAN(X/SQRT(1-X*X))
  216.       END; {ARCCOS}
  217.  
  218.    FUNCTION ARCCOSD (X:real) : real;
  219.       BEGIN
  220.          ARCCOSD := 57.295779513*ARCCOS(X)
  221.       END;
  222.  
  223.  
  224.    PROCEDURE Frame(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY: Integer);
  225.    VAR
  226.       i: Integer;
  227.    BEGIN
  228.       GotoXY(UpperLeftX, UpperLeftY);  Write('┌');
  229.       FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('─');
  230.       Write('┐');
  231.       FOR i:=UpperLeftY+1 TO LowerRightY-1 DO
  232.       BEGIN
  233.          GotX+1 TO LowerRightX-1 DO Write('─');
  234.       Write('┐');
  235.       FOR i:=UpperLeftY+1 TO LowerRightY-1 DO
  236.       BEGIN
  237.          GotoXY(UpperLeftX , i);  Write('│');
  238.          GotoXY(LowerRightX, i);  Write('│');
  239.       END;
  240.             GotoXY(UpperLeftX, LowerRightY);
  241.       Write('└');
  242.       FOR i:=UpperLeftX+1 TO LowerRightX-1 DO Write('─');
  243.       Write('┘');
  244.    END  { Frame };
  245.  
  246.  
  247. {      Astronomy Functions     }
  248.  
  249.  
  250.    FUNCTION Minutes (time :real) : integer;
  251.        BEGIN
  252.             Minutes := ABS(TRUNC(60*(FRAC(time))))
  253.        END;
  254.  
  255.    FUNCTION Seconds (time :real) : integer;
  256.        BEGIN
  257.             Seconds := ABS(TRUNC(60*(FRAC(60*FRAC(time)))))
  258.        END;
  259.  
  260.    PROCEDURE deg_min_sec (VAR angle:real);
  261.    BEGIN
  262.       IF NOT Display_Degrees_in_Decimal
  263.          THEN
  264.             BEGIN
  265.                IF angle < 0 THEN write('-') ELSE write(' ');
  266.                write(abs(trunc(angle)):3,'° ',minutes(angle):2,''' ',
  267.                      seconds(angle):2,'"');
  268.             END
  269.          ELSE
  270.             write(angle:8:4,'° ');
  271.        END;
  272.  
  273.    PROCEDURE hours_min_sec (VAR angle:real);
  274.    BEGIN
  275.       IF NOT Display_Degrees_in_Decimal
  276.          THEN
  277.             BEGIN
  278.                IF angle < 0 THEN write('-') ELSE write(' ');
  279.                write(abs(trunc(angle/15)):3,'h ',minutes(angle/15):2,'m ',
  280.                      seconds(angle/15):2,'s');
  281.             END
  282.          ELSE
  283.             write(angle/15:8:4,'h');
  284.        END;
  285.  
  286.    FUNCTION Daylight_Savings : boolean;
  287.      BEGIN
  288.        Daylight_Savings := (((Month + Day/100) < Daylight_savings_end)  AND
  289.                            ((Month + Day/100) > Daylight_savings_start) AND
  290.                            Daylight_savings_Correction_Enabled);
  291.      END;
  292.  
  293.    FUNCTION Anomaly (M,Ecc : real) : real;
  294.        VAR
  295.           d : real;
  296.           E : real;
  297.        BEGIN
  298.           E := M;
  299.           d := E - Ecc * SIN(E) - M;
  300.           WHILE abs(d) > 0.000001 DO
  301.              BEGIN
  302.                 E := E - d/(1-Ecc * COS(E));
  303.                 d := E - Ecc * SIN(E) - M
  304.              END;
  305.           Anomaly := 2 * ArcTan(sqrt((1+Ecc)/(1-Ecc)) * Tan(E/2))
  306.         END;{Anomaly}
  307.  
  308.    FUNCTION JulianDay (month,day,year :integer) : real;
  309.        VAR
  310.           A,B,C,D : real;
  311.  
  312.        BEGIN END;{Anomaly}
  313.  
  314.    FUNCTION JulianDay (month,day,year :integer) : real;
  315.        VAR
  316.           A,B,C,D : real;
  317.  
  318.        BEGIN
  319.          IF (month=1) OR (month=2)
  320.             THEN
  321.                BEGIN
  322.                   year := year-1;
  323.                   month := month+12
  324.                END;
  325.          IF (year>1582) OR (year=1582) AND (month>10) OR (year=1582) AND (month=10) AND (day>15)
  326.             THEN
  327.                BEGIN
  328.                   A := INT(year/100);
  329.                   B := 2-A+INT(A/4)
  330.                END
  331.             ELSE
  332.                BEGIN
  333.                   B := 0
  334.                END;
  335.  
  336.         C := INT(365.25*year);
  337.         D := INT(30.6001*(month+1));
  338.         JulianDay := B+C+D+day+1720994.5;
  339.       END;{JulianDay}
  340.  
  341. FUNCTION ECL_TO_RA (L,B : real) : real;
  342.    VAR
  343.      x : real;
  344.      y : real;
  345.      z : real;
  346.      zmax : real;
  347.    BEGIN
  348.      x := COSD(L);
  349.      y := SIND(L) * 0.91746406 - TAND(B) * 0.397818676;
  350.      z := ArcTanYX(y,x);
  351.      ECL_TO_RA := Z;
  352.    END;{ECL_TO_RA}
  353.  
  354. FUNCTION ECL_TO_DEC (L,B : real) : real;
  355.    BEGIN
  356.      ECL_TO_DEC := ARCSIND(SIND(B) * 0.91746406 + COSD(B) * SIND(L)
  357.                    * 0.397818676);
  358.    END;
  359.  
  360. FUNCTION LST_TO_LCT (LST,Long : real;
  361.                      Year, Day_of_year : integer;
  362.                      Zone_Corr : real) : real;
  363.    VAR
  364.      B,T,GST,GMT : real;
  365.    BEGIN
  366.         GST := LST + Long/15;
  367.         IF GST > 24 THEN GST := GST - 24;
  368.         IF GST < 0 THEN GST := GST + 24;
  369.         T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
  370.         B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
  371.              (24 * (year - 1900));
  372.         T := GST -  Day_of_Year  * 0.0657098 + B;
  373.         IF T > 24 THEN T := T - 24;
  374.         IF T < 0 THEN T := T + 24;
  375.         GMT := T * 0.99727;
  376.         LST_TO_LCT := GMT - Zone_Corr;
  377.     END;
  378.  
  379. PROCEDURE make_degrees_in_range(VAR n : real);
  380. BEGIN
  381.    WHILE n > 360 DO
  382.       n := n - 360;
  383.    WHILE n < 0 DO
  384.       n := n + 360;
  385. END;
  386.  
  387. PROCEDURE Make_Hours_in_Range(VAR n : real);
  388.    BEGIN
  389.       IF n > 24 THEN n := n - 24;
  390.       IF n <  0 THEN n := n + 24;
  391.    END;
  392.  
  393. PROCEDURE time_window;
  394.    BEGIN
  395.    WINDOW(1,1,80,25);
  396.    frame(56,1,80,10);
  397.    WINDOW(58,2,78,10);
  398.    GotoXY(1,1);
  399.    writeln('    ',month,'-',day,'-',year);
  400.    JD := Julianday(month,day,year);
  401.    daylight_savings_start :=
  402.           4 + (30 - Round(7 * Frac((Julianday(4,30,year) + 1.5)/7))) * 0.01;
  403.    daylight_savings_end :=
  404.           10 + (31 - Round(7 * Frac((Julianday(10,31,year) + 1.5)/7))) * 0.01;
  405.    writeln(' JD =',JD:10:0);
  406.    write('Long = ');
  407.    deg_min_sec(Longitude);
  408.    writeln;
  409.    write('Lat  = ');
  410.    deg_min_sec(Latitude);
  411.    writeln;
  412.    LCT := Hour + Minute/60 + 1/120;
  413.    writeln('    LCT = ',Hour:2,'h ',Minute:2,'m');
  414.    GMT := Zone_correction + LCT;
  415.    IF  Daylight_Savings THEN
  416.       GMT := GMT - 1;
  417.    Make_Hours_in_Range(GMT);
  418.    writeln ('    GMT = ',trunc(GMT):2,'h ',minutes(GMT):2,'m');
  419.    T := (JulianDay(1,0,year) - 2415020.0)/36525.0;
  420.    B := 24 - 6.6460656 - (2400.051262 * T) - (0.00002581 * T * T) +
  421.         (24 * (year - 1900));
  422.    GST := 0.0657098 * Day_of_year - B + GMT * 1.002738;
  423.    Make_Hours_in_Range(GST);
  424.    writeln ('    GST = ',trunc(GST):2,'h ',minutes(GST):2,'m');
  425.    LST := GST - Longitude/15;
  426.    Make_Hours_in_Range(LST);
  427.    writeln ('    LST = ',trunc(LST):2,'h ',minutes(LST):2,'m');
  428.    END;{time_window}
  429.  
  430. PROCEDURE Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n:integer);
  431.    BEGIN
  432.     WITH planets[n] DO
  433.      BEGIN
  434.        planet_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
  435.                                   - Julianday(1,0,Epoch)) / period;
  436.        Make_Degrees_in_Range(planet_long_of_ascend_node);
  437.        planet_mean_anomaly := planet_long_of_ascend_node + long_at_epoch
  438.                               - long_of_per;
  439.        planet_true_anomaly := DEGREE(ANOMALY(RADIAN(planet_mean_anomaly),Ecc));
  440.        planet_heliocentric_long := planet_true_anomaly + long_of_per;
  441.        Make_Degrees_in_Range(planet_heliocentric_long);
  442.        planet_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(planet_true_anomaly));
  443.        planet_helio_ecliptic_lat := ARCSIND(SIND(planet_heliocentric_long -
  444.                                     long_of_ascend_node) * SIND(Inc));
  445.     END;
  446.  
  447.    END;
  448.  
  449.  
  450. PROCEDURE Plot_planet(n,scale : integer);
  451. FUNCTION compute_radius_vector(n,angle : integer): real ;
  452. BEGIN
  453.    compute_radius_vector := planets[n].axis *
  454.                             (1 - planets[n].ecc * planets[n].ecc)
  455.                             /(1 + planets[n].ecc * COSD(angle
  456.                             - planets[n].long_of_per));
  457. END;
  458. VAR
  459. delX,delY,Xold,Xnew,Yold,Ynew,A : integer;
  460. BEGIN
  461.    Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(n);
  462.    Ynew := Round(scale * planet_radius_vector * SIND(planet_heliocentric_long))
  463.            + CenterY;
  464.    Xnew := Round(scale * Screen_Aspect * planet_radius_vector *
  465.            COSD(planet_heliocentric_long)) + CenterX;
  466.    FOR delX := -3 TO 3 DO
  467.       FOR delY := -1 TO 1 DO
  468.           plot(Xnew + delX,Ynew + delY,1);
  469.    GoToXY((Xnew DIV 8)+2,(Ynew DIV 8)+1);
  470.    write(planets[n].name);
  471.  
  472.    Xnew := CenterX + Round(compute_radius_vector(n,0) * scale * Screen_Aspect);
  473.    Ynew := CenterY;
  474.    FOR A := 1 TO 30 DO
  475.       BEGIN
  476.         + Round(compute_radius_vector(n,0) * scale * Screen_Aspect);
  477.    Ynew := CenterY;
  478.    FOR A := 1 TO 30 DO
  479.       BEGIN
  480.          Xold := Xnew;
  481.          Yold := Ynew;
  482.          Xnew := Round(compute_radius_vector(n,A*12) * scale * Screen_Aspect
  483.                  * COSD(A * 12)) + CenterX;
  484.          Ynew := Round(compute_radius_vector(n,A*12) * scale * SIND(A * 12))
  485.                  + CenterY;
  486.          Draw(Xold,Yold,Xnew,Ynew,1);
  487.       END;
  488. END;
  489.  
  490. {    Orbital Data for Planets   }
  491.  
  492.  
  493. PROCEDURE Orbital_Data;
  494. BEGIN
  495. WITH planets[1] DO
  496.    BEGIN
  497.       name :='Mercury   ';
  498.       Epoch := 1980;
  499.       Period := 0.24085;
  500.       Long_at_Epoch := 231.2973;
  501.       Long_of_Per := 77.1442128;
  502.       Ecc := 0.2056306;
  503.       Axis := 0.3870986;
  504.       Inc := 7.0043579;
  505.       Long_of_Ascend_Node := 48.0941733;
  506.       Size := 6.74;
  507.       Brightness := 0.000001918;
  508.    END;
  509. WITH planets[2] DO
  510.    BEGIN
  511.       name :='Venus     ';
  512.       Epoch := 1980;
  513.       Period := 0.61521;
  514.       Long_at_Epoch := 355.73352;
  515.       Long_of_Per := 131.2895792;
  516.       Ecc := 0.0067826;
  517.       Axis := 0.7233316;
  518.       Inc := 3.3944359;
  519.       Long_of_Ascend_Node := 76.4997524;
  520.       Size := 16.92;
  521.       Brightness := 0.00001721;
  522.    END;
  523. WITH planets[3] DO
  524.    BEGIN
  525.       name :='Earth     ';
  526.       Epoch := 1980;
  527.       Period := 1.00004;
  528.       Long_at_Epoch := 98.833540;
  529.       Long_of_Per := 102.596403;
  530.       Ecc := 0.016718;
  531.       Axis := 1.0;
  532.       Inc := 0;
  533.       Long_of_Ascend_Node := 0;
  534.       Size := 17;
  535.       Brightness := 0;
  536.    END;
  537. WITH planets[4] DO
  538.    BEGIN
  539.       name :='Mars      ';
  540.       Epoch := 1980;
  541.       Period := 1.88089;
  542.       Long_at_Epoch := 126.30783;
  543.       Long_of_Per := 335.6908166;
  544.       Ecc := 0.0933865;
  545.       Axis := 1.5236883;
  546.       Inc := 1.8498011;
  547.       Long_of_Ascend_Node := 49.4032001;
  548.       Size := 9.36;
  549.       Brightness := 0.00000454;
  550.    END;     Inc := 1.8498011;
  551.       Long_of_Ascend_Node := 49.4032001;
  552.       Size := 9.36;
  553.       Brightness := 0.00000454;
  554.    END;
  555. WITH planets[5] DO
  556.    BEGIN
  557.       name :='Jupiter   ';
  558.       Epoch := 1980;
  559.       Period := 11.86224;
  560.       Long_at_Epoch := 146.966365;
  561.       Long_of_Per := 14.0095493;
  562.       Ecc := 0.0484658;
  563.       Axis := 5.202561;
  564.       Inc := 1.3041819;
  565.       Long_of_Ascend_Node := 100.2520175;
  566.       Size := 196.74;
  567.       Brightness := 0.0001994;
  568.    END;
  569. WITH planets[6] DO
  570.    BEGIN
  571.       name :='Saturn    ';
  572.       Epoch := 1980;
  573.       Period := 29.4571;
  574.       Long_at_Epoch := 165.322242;
  575.       Long_of_Per := 92.6653974;
  576.       Ecc := 0.0556155;
  577.       Axis := 9.554747;
  578.       Inc := 2.4893741;
  579.       Long_of_Ascend_Node := 113.4888341;
  580.       Size := 165.60;
  581.       Brightness := 0.0001740;
  582.    END;
  583.  
  584. WITH planets[7] DO
  585.    BEGIN
  586.       name :='Uranus    ';
  587.       Epoch := 1980;
  588.       Period := 84.01247;
  589.       Long_at_Epoch := 228.0708551;
  590.       Long_of_Per := 172.7363288;
  591.       Ecc := 0.0463232;
  592.       Axis := 19.21814;
  593.       Inc := 0.7729895;
  594.       Long_of_Ascend_Node := 73.8768642;
  595.       Size := 65.80;
  596.       Brightness := 0.00007768;
  597.    END;
  598. WITH planets[8] DO
  599.    BEGIN
  600.       name :='Neptune   ';
  601.       Epoch := 1980;
  602.       Period := 164.79558;
  603.       Long_at_Epoch := 260.3578998;
  604.       Long_of_Per := 47.8672148;
  605.       Ecc := 0.0090021;
  606.       Axis := 30.10957;
  607.       Inc := 1.7716017;
  608.       Long_of_Ascend_Node := 131.5606494;
  609.       Size := 62.20;
  610.       Brightness := 0.00007597;
  611.    END;
  612.  
  613. WITH planets[9] DO
  614.    BEGIN
  615.       name :='Pluto     ';
  616.       Epoch := 1980;
  617.       Period := 250.9;
  618.       Long_at_Epoch := 209.439;
  619.       Long_of_Per := 222.972;
  620.       Ecc := 0.25387;
  621.       Axis := 39.78459;
  622.      := 250.9;
  623.       Long_at_Epoch := 209.439;
  624.       Long_of_Per := 222.972;
  625.       Ecc := 0.25387;
  626.       Axis := 39.78459;
  627.       Inc := 17.137;
  628.       Long_of_Ascend_Node := 109.941;
  629.       Size := 8.2;
  630.       Brightness := 0.000004073;
  631.    END;
  632. END; {Orbital_Data}
  633.  
  634. PROCEDURE Show_Menu;
  635. BEGIN
  636.    WINDOW(1,1,80,25);
  637.    CLRSCR;
  638.    Time_Window;
  639.    WINDOW(1,1,80,25);
  640.    GotoXY(1,1);
  641.  
  642.    FOR num := 1 TO Number_of_Planets DO
  643.       BEGIN
  644.          WITH planets[num] DO
  645.             writeln(num,'  ',name);
  646.       END;
  647.    writeln;
  648.    writeln('D  Change Date/Time');
  649.    writeln('I  Plot Inner Planets');
  650.    writeln('L  Change Long/Lat');
  651.    writeln('M  Menu');
  652.    writeln('O  Plot Outer Planets');
  653.    writeln('S  Set_up');
  654.    writeln('Q  Quit');
  655.    writeln;
  656.    writeln('Enter command:');
  657. END; {Show_menu}
  658.  
  659. PROCEDURE Set_Up;
  660. BEGIN
  661.    writeln('Display Degrees in Decimal Degrees ( Y/N )? ');
  662.    Read(KBD,Ch);
  663.    IF (Ch = 'Y') OR (Ch = 'y')
  664.       THEN Display_Degrees_in_Decimal := true;
  665.    IF (Ch = 'N') OR (Ch = 'n')
  666.       THEN Display_Degrees_in_Decimal := false;
  667.    writeln('Correct for daylight savings ( Y/N )?');
  668.    Read(KBD,Ch);
  669.    IF (Ch = 'Y') OR (Ch = 'y')
  670.       THEN Daylight_Savings_Correction_Enabled := true;
  671.    IF (Ch = 'N') OR (Ch = 'n')
  672.       THEN Daylight_Savings_Correction_Enabled := false;
  673.    writeln('Display Hour Angle in Degrees ( Y/N )?');
  674.    Read(KBD,Ch);
  675.    IF (Ch = 'Y') OR (Ch = 'y')
  676.       THEN Display_Hour_Angle_in_Degrees := true;
  677.    IF (Ch = 'N') OR (Ch = 'n')
  678.       ,Ch);
  679.    IF (Ch = 'Y') OR (Ch = 'y')
  680.       THEN Display_Hour_Angle_in_Degrees := true;
  681.    IF (Ch = 'N') OR (Ch = 'n')
  682.       THEN Display_Hour_Angle_in_Degrees := false;
  683.    writeln('Display RA in Degrees ( Y/N )?');
  684.    Read(KBD,Ch);
  685.    IF (Ch = 'Y') OR (Ch = 'y')
  686.       THEN Display_RA_in_Degrees := true;
  687.    IF (Ch = 'N') OR (Ch = 'n')
  688.       THEN Display_RA_in_Degrees := false;
  689. END;
  690.  
  691.  
  692. PROCEDURE Plot_Inner_Planets;
  693. BEGIN
  694.    HiRes;
  695.    HiResColor(1);
  696.    Window(1,1,80,25);
  697.    GotoXY(62,1);
  698.    writeln('    ',month,'-',day,'-',year);
  699.    GotoXY(62,2);
  700.    writeln(' JD =',JD:10:0);
  701.    plot(CenterX,CenterY,1);
  702.    FOR num := 1 TO 4 DO
  703.       plot_planet(num,scale_for_inner_planets);
  704. END;
  705.  
  706.  
  707. PROCEDURE Plot_Outer_Planets;
  708. BEGIN
  709.    HiRes;
  710.    HiResColor(1);
  711.    Window(1,1,80,25);
  712.    GotoXY(62,1);
  713.    writeln('    ',month,'-',day,'-',year);
  714.    GotoXY(62,2);
  715.    writeln(' JD =',JD:10:0);
  716.    plot(CenterX,CenterY,1);
  717.    plot_planet(3,scale_for_outer_planets);
  718.    FOR num := 5 TO 9 DO
  719.       plot_planet(num,scale_for_outer_planets);
  720. END;
  721.  
  722.  
  723. PROCEDURE Change_Date_Time;
  724.      BEGIN
  725.        time_window;
  726.        WINDOW(1,1,80,25);
  727.        writeln;
  728.        writeln ('Enter month  day  year');
  729.        readln (month,day,year);
  730.        writeln ('Enter hour minute');
  731.        readln (hour,minute);
  732.        day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
  733.        ur minute');
  734.        readln (hour,minute);
  735.        day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
  736.        time_window;
  737.      END;
  738.  
  739.  
  740.  
  741. PROCEDURE Change_Long_Lat;
  742.     BEGIN
  743.       time_window;
  744.       WINDOW(1,1,80,25);
  745.       writeln;
  746.       writeln ('Enter Longitude  Latitude');
  747.       readln (longitude,Latitude);
  748.       Zone_correction := trunc(longitude/15) + 1;
  749.       writeln ('Enter zone correction (Default = ',zone_correction,' ):');
  750.       readln (zone_correction);
  751.       time_window;
  752.     END;
  753.  
  754.  
  755. PROCEDURE Locate_Planet;
  756. LABEL
  757. End_of_Locate_Planet;
  758. BEGIN
  759.    Time_window;
  760.    WINDOW(1,1,80,25);
  761.    frame(1,1,54,9);
  762.    WINDOW(3,2,53,10);
  763.    GotoXY(1,1);
  764.    Locate_Position_of_Planet_in_Its_Own_Orbital_Plane(num);
  765.    WITH planets[num] DO
  766.      BEGIN
  767.        writeln('                          ',name);
  768.        write('long of ascend node  = ');
  769.        deg_min_sec(planet_long_of_ascend_node);
  770.        writeln;
  771.        write('mean anomaly         = ');
  772.        deg_min_sec(planet_mean_anomaly);
  773.        writeln;
  774.        write('true anomaly         = ');
  775.        deg_min_sec(planet_true_anomaly);
  776.        writeln;
  777.        write('heliocentric long    = ');
  778.        deg_min_sec(planet_heliocentric_long);
  779.        writeln;
  780.        writeln('radius vector        = ',planet_radius_vector:8:3,' AU');
  781.        write('helio ecliptic lat.  =');
  782.        deg_min_sec(planet_('radius vector        = ',planet_radius_vector:8:3,' AU');
  783.        write('helio ecliptic lat.  =');
  784.        deg_min_sec(planet_helio_ecliptic_lat);
  785.        writeln;
  786.      END;
  787.  
  788. {  Locate Position of Earth in Its Own Orbital Plane   }
  789.  
  790.    IF num = 3 THEN GoTO End_of_Locate_Planet;
  791.    WINDOW(40,2,53,10);
  792.    GotoXY(1,1);
  793.    writeln('     Earth');
  794.    WITH planets[3] DO
  795.      BEGIN
  796.  
  797.        earth_long_of_ascend_node := (360/365.2422) * (Julianday(month,day,year)
  798.                                   - Julianday(1,0,Epoch))/period;
  799.        Make_Degrees_in_Range(earth_long_of_ascend_node);
  800.        deg_min_sec(earth_long_of_ascend_node);
  801.        writeln;
  802.  
  803.        earth_mean_anomaly := earth_long_of_ascend_node + long_at_epoch
  804.                               - long_of_per;
  805.        deg_min_sec(earth_mean_anomaly);
  806.        writeln;
  807.  
  808.        earth_true_anomaly := DEGREE(ANOMALY(RADIAN(earth_mean_anomaly),Ecc));
  809.        deg_min_sec(earth_true_anomaly);
  810.        writeln;
  811.  
  812.        earth_heliocentric_long := earth_true_anomaly + long_of_per;
  813.        Make_Degrees_in_Range(earth_heliocentric_long);
  814.        deg_min_sec(earth_heliocentric_long);
  815.        writeln;
  816.  
  817.        earth_radius_vector := axis*(1-Ecc*Ecc)/(1+Ecc*COSD(earth_true_anomaly));
  818.        writeln(earth_radius_vector:8:3,' AU');
  819.        zero := 0;
  820.        deg_min_sec(zero);
  821.      END;
  822.  
  823.  
  824. {  Project Position of Planet onto plane of ecliptic   }
  825.  
  826.  
  827.    WINDOW(1,1,80,25);
  828.    frame(1,10,54,13);
  829.    WINDOW( 12,11,79,23);
  830.    GotoXY(1,1);
  831.    WITH planets[num] DO
  832.      BEGIN
  833.         y := SIND(planet_heliocentric_long - long_of_ascend_node) * COSD(Inc);
  834.         x := COSD(planet_heliocenO
  835.      BEGIN
  836.         y := SIND(planet_heliocentric_long - long_of_ascend_node) * COSD(Inc);
  837.         x := COSD(planet_heliocenO
  838.      BEGIN
  839.         y := SIND(planet_heliocentric_long - long_of_ascend_node) * COSD(Inc);
  840.         x := COSD(planet_heliocentric_long - long_of_ascend_node);
  841.         z :=ArcTanYX(y,x);
  842.         projected_heliocentric_long := Z + long_of_ascend_node;
  843.         Make_Degrees_in_Range(projected_heliocentric_long);
  844.         write('Projected Longitude  = ');
  845.         deg_min_sec(projected_heliocentric_long);
  846.         writeln;
  847.         Projected_Radius_vector := planet_radius_vector *
  848.                                    COSD(planet_helio_ecliptic_lat);
  849.         writeln('Projected Radius     = ',Projected_Radius_Vector:8:3,' AU');
  850.       END;
  851.  
  852. {    Calculate Ecliptical Coordinates   }
  853.  
  854.  
  855.    WINDOW(1,1,80,25);
  856.    frame(1,18,25,24);
  857.    WINDOW(3,19,24,23);
  858.    GotoXY(1,1);
  859.    IF (Planet_Radius_Vector < Earth_Radius_Vector)
  860.    THEN
  861.       Geocentric_Ecliptic_Long := Earth_Heliocentric_long + 180 +
  862.            ArcTanD(Projected_Radius_Vector * SIND(Earth_Heliocentric_long -
  863.            Projected_heliocentric_long)/(Earth_Radius_Vector -
  864.            Projected_Radius_Vector * COSD(Earth_Heliocentric_long -
  865.            Projected_heliocentric_long)))
  866.    ELSE
  867.       Geocentric_Ecliptic_Long := Projected_Heliocentric_Long +
  868.            ArcTanD(Earth_Radius_Vector * SIND(Projected_heliocentric_long -
  869.            Earth_Heliocentric_Long) / (Projected_Radius_Vector -
  870.            Earth_Radius_VeSIND(Projected_heliocentric_long -
  871.            Earth_Heliocentric_Long) / (Projected_Radius_Vector -
  872.            Earth_Radius_Vector * COSD(Projected_heliocentric_long -
  873.            Earth_Heliocentric_Long)));
  874.     Make_Degrees_in_Range(Geocentric_Ecliptic_Long);
  875.    writeln('      ECLIPTIC');
  876.    writeln;
  877.    write(' Long = ');
  878.    deg_min_sec(Geocentric_Ecliptic_Long);
  879.    writeln;
  880.    Geocentric_latitude := ArcTanD(Projected_Radius_Vector *
  881.                           TAND(planet_helio_ecliptic_lat) *
  882.                           SIND(Geocentric_Ecliptic_Long -
  883.                           Projected_Heliocentric_Long) / (Earth_Radius_Vector *
  884.                           SIND(Projected_Heliocentric_Long -
  885.                           Earth_Heliocentric_Long)));
  886.    write(' Lat  = ');
  887.    deg_min_sec(Geocentric_latitude);
  888.    writeln;
  889.  
  890.  
  891. {     Calculate Equatorial Coordinates    }
  892.  
  893.    WINDOW(1,1,80,25);
  894.    frame(27,18,52,24);
  895.    WINDOW(29,19,51,24);
  896.    GotoXY(1,1);
  897.    writeln('     EQUATORIAL');
  898.    writeln;INDOW(1,1,80,25);
  899.    frame(27,18,52,24);
  900.    WINDOW(29,19,51,24);
  901.    GotoXY(1,1);
  902.    writeln('     EQUATORIAL');
  903.    writeln;
  904.    RA := ECL_TO_RA(Geocentric_Ecliptic_Long,Geocentric_Latitude);
  905.    write(' RA  = ');
  906.    IF Display_RA_in_Degrees
  907.    THEN
  908.       deg_min_sec(RA)
  909.    ELSE
  910.       hours_min_sec(RA);
  911.    writeln;
  912.    DEC := ECL_TO_DEC(Geocentric_Ecliptic_Long,Geocentric_Latitude);
  913.    write(' DEC = ');
  914.    deg_min_sec(DEC);
  915.    writeln;
  916.    Hour_Angle := LST - RA/15;
  917.    IF Hour_Angle < 0 THEN Hour_Angle := Hour_Angle + 24;
  918.    Hour_Angle_in_Degrees := Hour_Angle * 15;
  919.    write(' HA  = ');
  920.    IF Display_Hour_Angle_in_Degrees
  921.    THEN
  922.       deg_min_sec(Hour_Angle_in_Degrees)
  923.    ELSE
  924.       hours_min_sec(Hour_Angle_in_Degrees);
  925.  
  926. {    Calculate Horizontal Coordinates   }
  927.  
  928.    WINDOW(1,1,80,25);
  929.    frame(54,18,80,24);
  930.    WINDOW(56,19,78,23);
  931.    GotoXY(1,1);
  932.    writeln('      HORIZONTAL');
  933.    writeln;
  934.    Altitude := ArcSinD(SIND(Dec) * SIND(Latitude) +
  935.                COSD(Dec) * COSD(Latitude) * COSD(Hour_Angle * 15));
  936.    write(' Alt  = ');
  937.    deg_min_sec(Altitude);
  938.    writeln;
  939.  
  940.    Azimuth := ArcCosD((SIND(Dec) - SIND(Latitude) * SIND(Altitude))/
  941.               (COSD(Latitude) * COSD(Altitude)));
  942.    IF SIND(Hour_Angle) > 0 THEN
  943.    Azimuth := 360 - Azimuth;
  944.    write(' Azim = ');
  945.    deg_min_sec(Azimuth);
  946.  
  947.  
  948. {   Calculate Time and Position of Rise and Set   }
  949.  
  950.  
  951.    WINDOW(1,1,80,25);
  952.    frame(1,14,54,17);
  953.    WINDOW(4,15,52,17);
  954.    GotoXY(1,1);
  955.    Azimuth_Rise := SIND(DEC)/COSD(Latitude);
  956.    IF (Azimuth_Rise < -1) OR (Azimuth_Rise > 1) THEN
  957.       IF Altitude < 0
  958.          THEN
  959.             writeln('never rises above horizon')
  960.          ELSE
  961.             writeln('never sets below horizon')
  962.    ELSE
  963.       BEGIN
  964.         Azimuth_Rise := ARCCOSD(Azimuth_Rise);
  965.         Azimuth_Set := 360 - Azimuth_Rise;
  966.         LST_Rise := 24 + RA/15 - (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
  967.         IF LST_Rise > 24 THEN LST_Rise := LST_Rise - 24;
  968.         LST_SET := RA/15 + (ARCCOSD(-TAND(Latitude) * TAND(DEC)))/15;
  969.         IF LST_SET > 24 THEN LST_SET := LST_SET - 24;
  970.  
  971.         LCT_SET := LST_TO_LCT(LST_SET,Longitude,Year,Day_of_Year,
  972.                     Zone_Correction);
  973.         IF Daylight_Savings THEN
  974.            LCT_SET := LCT_SET + 1;
  975.         IF LCT_SET < 0 THEN LCT_SET := LCT_SET + 24;
  976.  
  977.         LCT_RISE := LST_TO_LCT(LST_RISE,Longitude,Year,Day_of_Year,
  978.                     Zone_Correction);
  979.         IF Daylight_Savings THEN
  980.            LCT_RISE := LCT_RISE + 1;
  981.         IF LCT_RISE < 0 THEN LCT_RISE := LCT_RISE + 24;
  982.  
  983.         write('Rises at ',trunc(LCT_Rise):2,':',minutes(LCT_Rise):2,'  LCT');
  984.         write('      Azimuth  ');
  985.         deg_min_sec(Azimuth_Rise);
  986.         writeln;
  987.         write('Sets at  ',trunc(LCT_Set):2,':',minutes(LCT_Set):2,'  LCT');
  988.         write('      Azimuth  ');
  989.         deg_min_sec(Azimuth_Set);
  990.       END;
  991.  
  992.  
  993. {    Calculate Phase,Distance, Diameter, Magnitude   }
  994.  
  995.  
  996.    Window(1,1,80,25);
  997.    frame(56,11,80,17);
  998.    Window(57,12,79,16);
  999.    GotoXY(1,1);
  1000.  
  1001.    Phase := (1+COSD(Geocentric_ecliptic_long - planet_heliocentric_long))/2;
  1002.    writeln('Phase     = ',100 * Phase:6:2,'%');
  1003.    Distance_from_Earth := sqrt(Earth_Radius_Vector * Earth_Radius_Vector
  1004.                           + Planet_Radius_Vector * Planet_Radius_Vector
  1005.                           - 2 * Earth_Radius_Vector * Planet_Radius_Vector *
  1006.                           COSD(planet_heliocentric_long -
  1007.                           earth_heliocentric_long));
  1008.    writeln('Distance  = ',Distance_from_Earth:6:2,' AU');
  1009.  
  1010.    WITH Planets[num] DO
  1011.    Diameter := size / Distance_from_Earth;
  1012.    writeln('Diameter  = ',Diameter:6:2,'"');
  1013.  
  1014.    WITH Planets[num] DO
  1015.    Magnitude := 2.17147 * ln(Distance_from_Earth * Planet_Radius_Vector /
  1016.                 Brightness * sqrt(Phase)) - 26.7;
  1017.    writeln('Magnitude = ',Magnitude:6:2);
  1018. End_of_Locate_Planet:
  1019. END; {Locate_Planet}
  1020.  
  1021. BEGIN {Planets}
  1022.  
  1023. {      Initalize Default Parameters   }
  1024.  
  1025.    CLRSCR;
  1026.    MSDOS_time.AX := $2C00;{time call to MSDOS}
  1027.    MsDos(MSDOS_time);
  1028.    Hour :=Hi(MSDOS_time.CX);
  1029.    Minute := Lo(MSDOS_time.CX);
  1030.    MSDOS_time.AX := $2A00; {Date call to MSDOS}
  1031.    MsDos(MSDOS_time);
  1032.    Year := MSDOS_time.CX;
  1033.    Month := Hi(MSDOS_time.DX);
  1034.    Day := Lo(MSDOS_time.DX);
  1035.    Orbital_Data;
  1036.    day_of_year := TRUNC(Julianday(month,day,year)-Julianday(1,0,year));
  1037.  
  1038. {     Main Program     }
  1039.  
  1040.  
  1041. Show_Menu;
  1042. REPEAT
  1043.    Read(KBD, Ch);
  1044.    num := ord(Ch) - ord('0');
  1045.    TextMode;
  1046.    Bios_Display_Mode := $2D;
  1047.    ClrScr;
  1048.    CASE Ch OF
  1049.       
  1050.    Read(KBD, Ch);
  1051.    num := ord(Ch) - ord('0');
  1052.    TextMode;
  1053.    Bios_Display_Mode := $2D;
  1054.    ClrScr;
  1055.    CASE Ch OF
  1056.       'I','i' :  Plot_Inner_Planets;
  1057.       'O','o' :  Plot_Outer_Planets;
  1058.       'D','d' :  Change_Date_Time;
  1059.       'L','l' :  Change_Long_Lat;
  1060.       'M','m' :  Show_Menu;
  1061.       'S','s' :  Set_Up;
  1062.       '1','2','3','4','5','6','7','8','9'  : Locate_Planet;
  1063.    END; {case}
  1064.    Window(1,1,80,25);
  1065.    GotoXY(1,25);
  1066.    write('(Q)uit (M)enu (I)nner (O)uter (D)ate (L)ong   (1-9) Planet');
  1067. UNTIL (Ch = 'Q') OR (Ch = 'q');
  1068.    TextMode;
  1069.    Bios_Display_Mode := $2D;
  1070.    ClrScr;
  1071. END.
  1072. UNTIL (Ch = 'Q') OR (Ch = 'q');
  1073.    TextMode;
  1074.    Bios_Display_Mode := $2D;
  1075.    ClrScr;
  1076. END.